7-5 Days alive and free of invasive or non-invasive ventilation

This outcome was defined as the number of days alive and free of invasive or non-invasive ventilation from randomisation to 28 days. All patients dying within 28 days will be assigned zero free days.

Authors
Affiliations
James Totterdell

University of Sydney

Rob Mahar

University of Melbourne

Published

July 12, 2022

Preamble

Load packages
library(tidyverse)
library(labelled)
library(kableExtra)
library(cmdstanr)
library(posterior)
library(bayestestR)
library(bayesplot)
library(matrixStats)
library(plotly)
library(lubridate)
library(ggdist)
library(patchwork)

theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank()))
bayesplot_theme_set(theme_set(theme_classic(base_size = 10, base_family = "Palatino") +
  theme(panel.grid = element_blank(),
        strip.background = element_blank())))

color_scheme_set("red")
options(digits = 4)
Prepare dataset
devtools::load_all()
all_dat <- read_all_no_daily()
all_dat_daily <- read_all_daily()

acs_itt_dat <- all_dat %>% 
  filter_acs_itt() %>%
  transmute_model_cols_grp_aus_nz()
acs_itt_nona_dat <- acs_itt_dat %>% 
  filter(!is.na(out_dafv))

# Concurrent enrolments for C2
acs_itt_concurc2_dat <- acs_itt_dat %>%
  filter_concurrent_intermediate()
acs_itt_concurc2_nona_dat <- acs_itt_concurc2_dat %>% 
  filter(!is.na(out_dafv))

# Concurrent enrolments for C3
acs_itt_concurc3_dat <- acs_itt_dat %>%
  filter_concurrent_std_aspirin()
acs_itt_concurc3_nona_dat <- acs_itt_concurc3_dat %>% 
  filter(!is.na(out_dafv))

# Concurrent enrolments for C4
acs_itt_concurc4_dat <- acs_itt_dat %>%
  filter_concurrent_therapeutic()
acs_itt_concurc4_nona_dat <- acs_itt_concurc4_dat %>% 
  filter(!is.na(out_dafv))
Load models
ordmod0 <- cmdstan_model(
  "stan/ordinal/logistic_cumulative.stan") # No epoch or site
ordmod_epoch <- cmdstan_model(
  "stan/ordinal/logistic_cumulative_epoch.stan") # Epoch only
ordmod <- cmdstan_model(
  "stan/ordinal/logistic_cumulative_epoch_site.stan") # Full model
Code
make_summary_table <- function(dat, format = "html") {
  tdat <- dat %>%
    group_by(CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3", "C4"),
      labels = intervention_labels2()$CAssignment[-1])) %>%
    summarise(
      Patients = n(),
      Known = sum(!is.na(out_dafv)),
      Deaths = sprintf("%i (%.0f%%)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      `Any ventilation` = sprintf("%i (%.0f%%)", 
                                  sum(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE),
                                  100 * mean(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE)),
      `DAFV, Median (Q1, Q3)` = sprintf(
        "%.0f (%.0f, %.0f)", 
        median(out_dafv, na.rm = T), 
        quantile(out_dafv, 0.25, na.rm = TRUE), 
        quantile(out_dafv, 0.75, na.rm = TRUE))
    ) %>%
    bind_rows(
      dat %>%
    group_by(CAssignment = "Overall") %>%
    summarise(
      Patients = n(),
      Known = sum(!is.na(out_dafv)),
      Deaths = sprintf("%i (%.0f%%)", sum(D28_death, na.rm = TRUE), 100 * mean(D28_death, na.rm = TRUE)),
      `Any ventilation` = sprintf("%i (%.0f%%)", 
                              sum(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE),
                              100 * mean(D28_OutcomeDaysFreeOfVentilation < 28, na.rm = TRUE)),
      `DAFV, Median (Q1, Q3)` = sprintf(
        "%.0f (%.0f, %.0f)", 
        median(out_dafv, na.rm = T), 
        quantile(out_dafv, 0.25, na.rm = TRUE), 
        quantile(out_dafv, 0.75, na.rm = TRUE))
    )
    ) %>%
    rename(`Anticoagulation\nintervention` = CAssignment)
  kable(
    tdat,
    format = format,
    align = "lrrrrr",
    booktabs = TRUE,
    linesep = ""
  ) %>%
    kable_styling(
      font_size = 9,
      latex_options = "HOLD_position"
    ) %>%
    row_spec(nrow(tdat), bold = T)
}


make_primary_model_data <- function(
    dat, 
    vars = NULL,
    beta_sd_var = NULL,
    ctr = contr.orthonorm,
    p_mult = 2) {
  
  X <- make_X_design(dat, vars, ctr)
  attX <- attr(X, "contrasts")
  X <- X[, -1]
  attr(X, "contrasts") <- attX
  nXtrt <- sum(grepl("rand", colnames(X)))
  
  beta_sd <- c(rep(1, nXtrt), beta_sd_var)
  
  epoch <- dat$epoch
  M_epoch <- max(dat$epoch)
  region <- dat$ctry_num
  M_region <- max(region)
  site <- dat$site_num
  M_site <- max(site)
  region_by_site <- dat %>% 
    dplyr::count(ctry_num, site_num) %>% 
    pull(ctry_num)

  y_raw <- ordered(dat$out_dafv)
  y_mod <- as.integer(y_raw)
  N <- dim(X)[1]
  K <- dim(X)[2]
  unique_y <- length(unique(y_mod))
  list(N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
       M_region = M_region,
       M_site = M_site, site = site,
       M_epoch = M_epoch, epoch = epoch,
       region_by_site = region_by_site,
       p_par = p_mult * rep(1 / unique_y, unique_y),
       beta_sd = beta_sd)  
}

Outcome Derivation

This outcome was defined as the number of days alive and free of positive pressure ventilation (invasive or non-invasive) to 28 days. All patients dying within 28 days will be assigned zero free days.

The outcome is calculated for a patient as:

  • missing, if day 28 mortality is unknown (D28_PatientStatus) or if the patient was known to have been alive at day 28 but number of days requiring ventilation (D28_OutcomeDaysFreeOfVentilation) is unknown
  • 0 if they died by day 28 (D28_PatientStatus)
  • min(28, D28_OutcomeDaysFreeOfVentilation) otherwise

Descriptive Analyses

Overall

The overall distribution of days alive and free of ventilation (DAFV) are reported in Figure 1 for all participants in the ACS-ITT set.

There were few participants who required any ventilation during hospital who did not die (DAFV of 0),

Code
pdat <- acs_itt_dat %>%
  dplyr::count(dafv = addNA(ordered(out_dafv, levels = 0:28)), .drop = F) %>%
  mutate(p = n / sum(n))
ggplot(pdat, aes(dafv, p)) +
  geom_col() +
  labs(
    x = "Days alive and free of ventilation", 
    y = "Proportion"
  )

Code
save_tex_table(
  make_summary_table(acs_itt_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-summary"))
make_summary_table(acs_itt_dat)
Anticoagulation intervention Patients Known Deaths Any ventilation DAFV, Median (Q1, Q3)
Low-dose 610 596 19 (3%) 34 (6%) 28 (28, 28)
Intermediate-dose 613 603 15 (2%) 23 (4%) 28 (28, 28)
Low-dose with aspirin 283 281 10 (4%) 18 (6%) 28 (28, 28)
Therapeutic-dose 50 50 6 (12%) 7 (14%) 28 (28, 28)
Overall 1556 1530 50 (3%) 82 (5%) 28 (28, 28)
Summaty of DAFV.

Anti-coagulation

Code
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    CAssignment = factor(CAssignment, labels = str_replace(intervention_labels()$CAssignment[-1], "<br>", "\n")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive.pdf"), p, height = 3, width = 6)
p

Age

DAFV by age
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    agegrp = cut(AgeAtEntry, c(18, seq(25, 75, 5), 100), include.lowest = T, right = F),
    dafv = fct_rev(ordered(out_dafv, levels = 0:28)), 
    .drop = F) %>%
  group_by(agegrp) %>%
  mutate(p = n / sum(n))
pdat2 <- pdat %>%
  group_by(agegrp) %>%
  summarise(n = sum(n))
p1 <- ggplot(pdat2, aes(agegrp, n)) +
  geom_col(colour = "grey40", fill = "grey40") +
  geom_vline(xintercept = 60, linetype = 2) +
  labs(y = "Number of\nparticipants",
       x = "Age at entry") +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45))
p2 <- ggplot(pdat, aes(agegrp, p)) +
  geom_col(aes(fill = dafv)) +
  labs(x = "Age", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  geom_vline(xintercept = 8.5, linetype = 2) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-age.pdf"), p, height = 2.5, width = 6)
p

Country

DAFV by country
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    country,
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(country) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
pdat2 <- pdat %>%
  group_by(country) %>%
  summarise(n = sum(n))

p1 <- ggplot(pdat2, aes(country, n)) +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Country of enrolment")

p2 <- ggplot(pdat, aes(country, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-country.pdf"), p, height = 2.5, width = 6)
p

Site

DAFV by site
pdat <- all_dat %>%
  filter_acs_itt() %>%
  filter(!is.na(out_dafv)) %>%
  dplyr::count(
    Country = factor(PT_CountryName, levels = c("India", "Australia", "Nepal", "New Zealand"),
                     labels = c("India", "Australia", "Nepal", "New\nZealand")),
    Site = fct_infreq(Location),
    dafv = ordered(out_dafv, levels = 0:28)) %>%
  complete(dafv = ordered(0:28), nesting(Country, Site), fill = list(n = 0)) %>%
  group_by(Country, Site) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup() %>%
  mutate(
    Country = droplevels(Country),
    Site = droplevels(Site)
  )
pdat2 <- pdat %>%
  group_by(Country, Site) %>%
  summarise(n = sum(n)) %>%
  ungroup()
p1 <- ggplot(pdat2, aes(Site, n)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        panel.border = element_rect(fill = NA))
p2 <- ggplot(pdat, aes(Site, p)) +
  facet_grid( ~ Country, scales = "free_x", space = "free_x") +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 1)) +
  theme(axis.text.x = element_text(hjust = 1, angle = 45)) +
  theme(legend.key.size = unit(0.25, "lines"))
p <- p1 / p2 +
  plot_layout(guides = 'collect')
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-country-site.pdf"), p, height = 4, width = 6.25)
p

Calendar Time

DAFV by calendar date
pdat <- acs_itt_nona_dat %>%
  dplyr::count(
    yr = year(RandDate), mth = month(RandDate),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(yr, mth) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p1 <- pdat %>%
  group_by(yr, mth) %>%
  summarise(n = sum(n)) %>%
  ggplot(., aes(mth, n))  +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
    geom_col() +
    labs(
      y = "Number of\nparticipants", 
      x = "Calendar date (month of year)") +
  scale_x_continuous(breaks = 1:12)
p2 <- ggplot(pdat, aes(mth, p)) +
  facet_grid( ~ yr, drop = T, scales = "free_x", space = "free") +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Calendar date (month of year)", y = "Cumulative\nproportion") +
  scale_fill_viridis_d("DAFV", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE)) +
  theme(legend.key.size = unit(0.25, "lines")) +
  scale_x_continuous(breaks = 1:12)
p <- p1 | p2
path <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(path, "7-5-calendar-time.pdf"), p, height = 2, width = 6)
p

Sample Cumulative Logits

Plot sample cumulative logits
trt_counts <- acs_itt_nona_dat %>%
  dplyr::count(CAssignment, out_dafv) %>%
  complete(CAssignment, out_dafv, fill = list(n = 0)) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n))
trt_logit <- trt_counts %>% 
  group_by(CAssignment) %>% 
  mutate(clogit = logit(cumsum(p))) %>%
  group_by(out_dafv) %>%
  mutate(rel_clogit = clogit - mean(clogit)) %>%
  filter(out_dafv != 28)

ggplot(trt_logit, aes(CAssignment, rel_clogit)) +
  facet_wrap( ~ out_dafv) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Plot sample cumulative logits
ggplot(trt_logit, aes(out_dafv, rel_clogit)) +
  facet_wrap( ~ CAssignment) +
  geom_point() +
  labs(y = "Relative (to mean) sample cumulative logit")

Modelling

The primary analysis uses the ACS-ITT set, so all participants have been randomised to the anticoagulation domain and some participants may have been randomised to the antiviral domain.

Primary Model

Fit primary model
fit_primary_model <- function(dat, ctr = contr.orthonorm, save = FALSE) {
  mdat <- make_primary_model_data(
    dat, c("inelgc3", "agegte60", "ctry"), 
    c(10, 2.5, 1, 1), 
    ctr)
  snk <- capture.output(
  mfit <- ordmod$sample(
    data = mdat,
    seed = 3510392,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0,
    adapt_delta = 0.98
  )) 
  if (save) {
    mfit$save_object(file.path("outputs", "models", "secondary", "7-5.rds"))    
  }
  mdrws <- as_draws_rvars(mfit$draws())
  # Add names
  epoch_map <- dat %>% dplyr::count(epoch, epoch_lab)
  site_map <- dat %>% dplyr::count(site_num, site)
  names(mdrws$beta) <- colnames(mdat$X)
  names(mdrws$gamma_epoch) <- epoch_map$epoch_lab
  names(mdrws$gamma_site) <- site_map$site
  # Convert to treatment log odds ratios
  mdrws$Acon <- attr(mdat$X, "contrasts")$randA %**% mdrws$beta[grepl("randA[0-9]", names(mdrws$beta))]
  mdrws$Ccon <- attr(mdat$X, "contrasts")$randC %**% mdrws$beta[grepl("randC[0-9]", names(mdrws$beta))]
  mdrws$trtA <- mdrws$Acon[-1] - mdrws$Acon[1]
  mdrws$trtC <- mdrws$Ccon[-1] - mdrws$Ccon[1]
  mdrws$compare <- c(
    "Intermediate vs low" = exp(mdrws$trtC[1]),
    "Low with aspirin vs low" = exp(mdrws$trtC[2]),
    "Therapeutic vs low" = exp(mdrws$trtC[3]),
    "Intermediate vs low with aspirin" = exp(mdrws$trtC[1] - mdrws$trtC[2]),
    "Intermediate vs therapeutic" = exp(mdrws$trtC[1] - mdrws$trtC[3]),
    "Low with aspirin vs therapeutic" = exp(mdrws$trtC[2] - mdrws$trtC[3])
  )
  mdrws$OR <- c(
    setNames(exp(mdrws$trtC), c("Intermediate", "Low with aspirin", "Therapeutic")),
    "Ineligible aspirin" = exp(mdrws$beta[which(grepl("inelg", colnames(mdat$X)))]),
    "Age \u2265 60" = exp(mdrws$beta[which(grepl("age", colnames(mdat$X)))]),
    setNames(exp(mdrws$beta[which(grepl("ctry", colnames(mdat$X)))]), c("AU/NZ", "NP"))
  )
  return(list(dat = mdat, fit = mfit, drws = mdrws))
}
res <- fit_primary_model(acs_itt_nona_dat, save = TRUE)
Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(res$drws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-summary-table")
odds_ratio_summary_table_rev(res$drws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.40 (0.81, 2.45) 1.46 (0.42) 0.88
Low with aspirin 1.16 (0.62, 2.19) 1.22 (0.41) 0.68
Therapeutic 0.41 (0.15, 1.15) 0.47 (0.27) 0.05
Ineligible aspirin 0.29 (0.09, 1.08) 0.37 (0.28) 0.03
Age ≥ 60 0.50 (0.31, 0.83) 0.52 (0.13) 0.00
AU/NZ 1.10 (0.27, 4.92) 1.47 (1.28) 0.55
NP 0.56 (0.15, 2.32) 0.74 (0.68) 0.20
Posterior summaries for model parameters (fixed-effects).
Odds ratio summary for epoch and site
p <- plot_epoch_site_terms(
  exp(res$drws$gamma_epoch),
  exp(res$drws$gamma_site),
  factor(res$dat$region_by_site, 
         labels = c("India", "Australia\nNew Zealand", "Nepal"))
)
pth <- file.path("outputs", "figures", "outcomes", "secondary", "7-5-primary-model-epoch-site-terms.pdf")
ggsave(pth, p, width = 6, height = 4.5)
p

Code
p <- plot_or_densities(res$drws$compare)
p

Code
res$fit$summary(c("alpha", "beta", "gamma_site", "gamma_epoch", "tau_site", "tau_epoch")) %>%
  print(n = Inf)
# A tibble: 68 × 10
   variable            mean   median    sd   mad      q5     q95  rhat ess_bulk
   <chr>              <dbl>    <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>
 1 alpha[1]        -4.37    -4.36    0.704 0.682 -5.55   -3.23    1.00    2883.
 2 alpha[2]        -4.32    -4.32    0.704 0.684 -5.50   -3.18    1.00    2869.
 3 alpha[3]        -4.30    -4.29    0.704 0.682 -5.48   -3.16    1.00    2858.
 4 alpha[4]        -4.28    -4.27    0.704 0.682 -5.46   -3.14    1.00    2850.
 5 alpha[5]        -4.26    -4.25    0.703 0.679 -5.43   -3.11    1.00    2842.
 6 alpha[6]        -4.23    -4.23    0.703 0.679 -5.41   -3.10    1.00    2840.
 7 alpha[7]        -4.19    -4.19    0.702 0.678 -5.37   -3.06    1.00    2832.
 8 alpha[8]        -4.17    -4.17    0.702 0.677 -5.35   -3.03    1.00    2823.
 9 alpha[9]        -4.13    -4.13    0.702 0.679 -5.30   -3.00    1.00    2812.
10 alpha[10]       -4.10    -4.09    0.701 0.677 -5.27   -2.97    1.00    2809.
11 alpha[11]       -4.08    -4.07    0.701 0.677 -5.24   -2.95    1.00    2805.
12 alpha[12]       -4.06    -4.05    0.701 0.678 -5.22   -2.93    1.00    2800.
13 alpha[13]       -4.02    -4.01    0.701 0.678 -5.19   -2.89    1.00    2793.
14 alpha[14]       -3.99    -3.98    0.700 0.677 -5.16   -2.86    1.00    2782.
15 alpha[15]       -3.94    -3.93    0.700 0.676 -5.11   -2.81    1.00    2778.
16 alpha[16]       -3.86    -3.85    0.699 0.674 -5.03   -2.74    1.00    2780.
17 alpha[17]       -3.81    -3.80    0.699 0.673 -4.97   -2.68    1.00    2771.
18 beta[1]         -0.753   -0.758   0.420 0.426 -1.44   -0.0600  1.00   18460.
19 beta[2]          0.118    0.117   0.226 0.228 -0.251   0.493   1.00   20145.
20 beta[3]          0.555    0.553   0.279 0.281  0.0984  1.02    1.00   18259.
21 beta[4]         -0.121   -0.115   0.694 0.688 -1.28    1.00    1.00   16789.
22 beta[5]          0.893    0.882   0.530 0.526  0.0392  1.79    1.00   23035.
23 beta[6]         -1.21    -1.23    0.633 0.630 -2.22   -0.134   1.00   24360.
24 beta[7]         -0.691   -0.692   0.255 0.254 -1.10   -0.269   1.00   19969.
25 beta[8]          0.112    0.0971  0.732 0.729 -1.07    1.35    1.00   10389.
26 beta[9]         -0.560   -0.578   0.692 0.676 -1.67    0.594   1.00    6683.
27 gamma_site[1]    2.07     1.82    1.50  1.27   0.139   4.83    1.00    8903.
28 gamma_site[2]    2.25     2.00    1.51  1.27   0.331   4.96    1.00    7889.
29 gamma_site[3]    1.48     1.25    1.54  1.31  -0.534   4.26    1.00   10176.
30 gamma_site[4]   -1.06    -1.08    0.601 0.572 -2.02   -0.0458  1.00    4793.
31 gamma_site[5]   -0.0208  -0.0425  0.596 0.565 -0.956   0.996   1.00    4206.
32 gamma_site[6]   -1.61    -1.62    0.626 0.610 -2.63   -0.582   1.00    5828.
33 gamma_site[7]    2.27     2.02    1.52  1.27   0.346   5.06    1.00    8546.
34 gamma_site[8]   -0.440   -0.450   0.692 0.666 -1.56    0.707   1.00    6550.
35 gamma_site[9]    1.47     1.40    0.848 0.808  0.221   2.97    1.00    7327.
36 gamma_site[10]   1.98     1.74    1.50  1.28   0.0214  4.74    1.00    8735.
37 gamma_site[11]  -0.387   -0.400   0.657 0.628 -1.43    0.716   1.00    5897.
38 gamma_site[12]  -0.569   -0.435   1.03  0.890 -2.41    0.933   1.00   16270.
39 gamma_site[13]   0.627    0.406   1.18  0.866 -0.898   2.78    1.00   19907.
40 gamma_site[14]  -1.05    -0.906   1.04  1.05  -2.94    0.361   1.00   11400.
41 gamma_site[15]   0.534    0.333   1.21  0.884 -1.09    2.70    1.00   19995.
42 gamma_site[16]  -0.989   -0.875   0.988 1.02  -2.78    0.369   1.00   11863.
43 gamma_site[17]  -0.976   -0.894   0.887 0.905 -2.57    0.253   1.00   11363.
44 gamma_site[18]   0.425    0.232   1.21  0.900 -1.23    2.61    1.00   21990.
45 gamma_site[19]   0.336    0.166   1.22  0.906 -1.39    2.49    1.00   21820.
46 gamma_site[20]  -0.587   -0.459   1.03  0.892 -2.44    0.902   1.00   16073.
47 gamma_site[21]   0.438    0.246   1.22  0.901 -1.24    2.63    1.00   22541.
48 gamma_site[22]   0.502    0.294   1.23  0.903 -1.15    2.71    1.00   21711.
49 gamma_site[23]   1.15     0.938   1.15  1.00  -0.233   3.29    1.00   14574.
50 gamma_site[24]   0.625    0.409   1.18  0.872 -0.924   2.80    1.00   18145.
51 gamma_site[25]  -0.270   -0.126   0.579 0.359 -1.39    0.454   1.00   11367.
52 gamma_site[26]  -0.00143  0.00867 0.557 0.335 -0.965   0.859   1.00   15002.
53 gamma_epoch[1]   0        0       0     0      0       0      NA         NA 
54 gamma_epoch[2]  -0.229   -0.169   0.402 0.332 -0.969   0.338   1.00   10447.
55 gamma_epoch[3]  -0.333   -0.268   0.478 0.427 -1.20    0.345   1.00    8933.
56 gamma_epoch[4]  -0.266   -0.224   0.464 0.415 -1.08    0.445   1.00    9924.
57 gamma_epoch[5]  -0.545   -0.474   0.561 0.551 -1.58    0.231   1.00    7721.
58 gamma_epoch[6]  -0.373   -0.325   0.554 0.507 -1.33    0.469   1.00   10405.
59 gamma_epoch[7]  -0.480   -0.430   0.569 0.550 -1.48    0.357   1.00    8708.
60 gamma_epoch[8]  -0.567   -0.515   0.577 0.564 -1.58    0.283   1.00    8024.
61 gamma_epoch[9]  -0.795   -0.751   0.608 0.617 -1.86    0.0821  1.00    6530.
62 gamma_epoch[10] -1.08    -1.05    0.680 0.698 -2.26   -0.0280  1.00    5706.
63 gamma_epoch[11] -0.595   -0.550   0.616 0.616 -1.66    0.328   1.00    7632.
64 gamma_epoch[12] -0.307   -0.281   0.645 0.586 -1.40    0.741   1.00   10145.
65 tau_site[1]      1.83     1.68    0.733 0.603  0.937   3.22    1.00    5368.
66 tau_site[2]      1.20     1.12    0.680 0.616  0.231   2.42    1.00    6674.
67 tau_site[3]      0.647    0.483   0.619 0.454  0.0438  1.80    1.00   11241.
68 tau_epoch        0.439    0.412   0.237 0.224  0.0916  0.863   1.00    5973.
# … with 1 more variable: ess_tail <dbl>
Code
res$fit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.7775 0.7944 0.7418 0.8248 0.7489 0.7897 0.8120 0.7979
Code
mcmc_trace(res$drws["beta"])

Code
mcmc_trace(res$drws["alpha"])

Code
mcmc_trace(res$drws["gamma_site"])

Code
mcmc_trace(res$drws["gamma_epoch"])

Posterior Predictive

Code
y_raw <- as.integer(levels(res$dat$y_raw))
y_ppc <- res$drws$y_ppc
y_ppc_raw <- rfun(\(x) y_raw[x])(y_ppc)
ppc_dat <- bind_cols(acs_itt_nona_dat, tibble(y_ppc = y_ppc_raw)) %>%
  mutate(CAssignment = factor(
    CAssignment, 
    levels = names(intervention_labels_short_break()$CAssignment),
    labels = intervention_labels_short_break()$CAssignment))

grp_ppc <- function(grp, d = 0:28) {
  ppc_dat %>%
    group_by(grp = {{grp}}) %>%
    summarise(
      y_lteq = map_dbl(d, ~ mean(out_dafv <= .x)),
      ypp_lteq = map(d, ~ rvar_mean(y_ppc <= .x)),
      y_eq = map_dbl(d, ~ mean(out_dafv == .x)),
      ypp_eq = map(d, ~ rvar_mean(y_ppc == .x))
    ) %>%
    unnest(c(ypp_lteq, ypp_eq)) %>%
    mutate(days = d) %>%
    ungroup() %>%
    pivot_longer(
      y_lteq:ypp_eq, 
      names_to = c("response", "event"),
      names_sep = "_",
      values_to = "posterior")
}

plot_grp_ppc <- function(dat, lab = "", xlab = "Probability") {
  ggplot(dat %>% 
           filter(response == "ypp", event == "eq"), 
         aes(y = days)) +
    facet_wrap( ~ grp, nrow = 1, scales = "free_x") +
    stat_slabinterval(aes(xdist = posterior), fatten_point = 1)  +
    geom_point(data = dat %>% filter(response == "y", event == "eq"), 
               aes(y = days, x = mean(posterior)),
               colour = "red",
               shape = 23) +
    scale_x_continuous(
      xlab, breaks = c(0, 0.3, 0.6),
      sec.axis = sec_axis(
        ~ . , name = lab, breaks = NULL, labels = NULL)) +
    scale_y_continuous("Days") +
    theme(strip.text = element_text(size = rel(0.7)),
          axis.title.x = element_text(size = rel(0.7)),
          axis.text.x = element_text(size = rel(0.65)),
          axis.title.y = element_text(size = rel(0.75)),
          axis.title.x.bottom = element_blank())
}

pp_C <- grp_ppc(CAssignment)
pp_ctry <- grp_ppc(country)
pp_epoch <- grp_ppc(epoch)
pp_site <- grp_ppc(site) %>%
  left_join(ppc_dat %>% dplyr::count(site, country), 
            by = c("grp" = "site"))

p1 <- plot_grp_ppc(pp_C, "Anticoagulation", "")
p2 <- plot_grp_ppc(pp_ctry, "Country", "") 
p3 <- plot_grp_ppc(pp_epoch, "Epoch", "")
p4 <- plot_grp_ppc(
  pp_site %>% filter(country == "IN"), "Sites India", "")
p5 <- plot_grp_ppc(
  pp_site %>% filter(country == "AU"), "Sites Australia", "")
p6 <- plot_grp_ppc(
  pp_site %>% filter(country == "NP"), "Sites Nepal", "")
p7 <- plot_grp_ppc(
  pp_site %>% filter(country == "NP"), "Sites New Zealand", "")

g1 <- (p1 | p2) / p3
g2 <- p4 / p5 / (p6 | p7)

pth1 <- file.path(
  "outputs", "figures", "outcomes", "secondary",
  "7-5-primary-model-acs-itt-ppc1.pdf")
pth2 <- file.path(
  "outputs", "figures", "outcomes", "secondary",
  "7-5-primary-model-acs-itt-ppc2.pdf")

ggsave(pth1, g1, width = 6, height = 5, device = cairo_pdf)
ggsave(pth2, g2, width = 6, height = 6, device = cairo_pdf)

g1

Code
g2

Sensitivity: Concurrent controls

Three separate analyses are conducted based on concurrent randomisations as was done for the primary outcome. For these models, the epoch term has been removed.

Intermediate dose
Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc2_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-concurrent-intermediate-summary"))
make_summary_table(acs_itt_concurc2_dat)
Anticoagulation intervention Patients Known Deaths Any ventilation DAFV, Median (Q1, Q3)
Low-dose 610 596 19 (3%) 34 (6%) 28 (28, 28)
Intermediate-dose 613 603 15 (2%) 23 (4%) 28 (28, 28)
Overall 1223 1199 34 (3%) 57 (5%) 28 (28, 28)
Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc2_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc2_nona_dat$out_dafv)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 2.5, 1, 1)
)
snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 75136,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(2) %**% mdrws$beta[1])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt), 
              setNames(exp(mdrws$beta[2:4]), c("Age \u2265 60", "Australia/New Zealand", "Nepal")))
Sample distribution
pdat <- acs_itt_concurc2_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2"),
      labels = c("Low", "Intermediate")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv)), colour = NA) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive-concurrent-intermediate.pdf"), p, height = 2.5, width = 4)
p

Posterior contrast
p <- plot_or_densities(mdrws$compare) +
  geom_text(aes(x = 1, y = 1 - 0.1, 
                label = sprintf("Pr(OR > 1) = %.2f", Pr(RV > 1))), 
            hjust = 0, nudge_x = 0.01, family = "Palatino") +
  scale_x_log10(breaks = c(0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6)) +
  labs(y = "Comparison")
p

Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-concurrent-intermediate-summary-table")
odds_ratio_summary_table_rev(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.51 (0.91, 2.59) 1.57 (0.43) 0.94
Age ≥ 60 0.43 (0.25, 0.74) 0.44 (0.13) 0.00
Australia/New Zealand 0.73 (0.36, 1.62) 0.80 (0.33) 0.21
Nepal 0.92 (0.40, 2.44) 1.05 (0.54) 0.43
Posterior summaries for model parameters.
Code
mfit$summary()
# A tibble: 40 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -3.78e+2 -3.78e+2 3.29e+0 3.22e+0 -3.84e+2 -3.73e+2  1.00    8006.
 2 beta_raw[…  2.94e-1  2.92e-1 1.90e-1 1.88e-1 -1.79e-2  6.09e-1  1.00   23243.
 3 beta_raw[… -3.40e-1 -3.42e-1 1.11e-1 1.11e-1 -5.23e-1 -1.58e-1  1.00   15990.
 4 beta_raw[… -3.00e-1 -3.12e-1 3.84e-1 3.84e-1 -9.12e-1  3.55e-1  1.00   20207.
 5 beta_raw[… -6.35e-2 -8.54e-2 4.57e-1 4.47e-1 -7.87e-1  7.28e-1  1.00   21016.
 6 p[1]        1.95e-2  1.92e-2 4.30e-3 4.25e-3  1.31e-2  2.71e-2  1.00   16213.
 7 p[2]        6.52e-4  4.69e-4 6.23e-4 4.59e-4  4.28e-5  1.89e-3  1.00   16237.
 8 p[3]        6.58e-4  4.65e-4 6.39e-4 4.57e-4  4.56e-5  1.94e-3  1.00   18456.
 9 p[4]        6.49e-4  4.63e-4 6.25e-4 4.48e-4  4.29e-5  1.89e-3  1.00   14798.
10 p[5]        6.56e-4  4.75e-4 6.24e-4 4.57e-4  4.78e-5  1.87e-3  1.00   17414.
# … with 30 more rows, and 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9209 1.0093 1.0153 0.9283 0.9173 1.0053 0.9836 0.9998
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])

Low-dose with aspirin

Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc3_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-concurrent-stdaspirin-summary"))
make_summary_table(acs_itt_concurc3_dat)
Anticoagulation intervention Patients Known Deaths Any ventilation DAFV, Median (Q1, Q3)
Low-dose 299 291 11 (4%) 20 (7%) 28 (28, 28)
Intermediate-dose 298 293 12 (4%) 15 (5%) 28 (28, 28)
Low-dose with aspirin 283 281 10 (4%) 18 (6%) 28 (28, 28)
Overall 880 865 33 (4%) 53 (6%) 28 (28, 28)
Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc3_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc3_nona_dat$out_dafv)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1)
)
snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 135356,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Low with aspirin vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs low with aspirin" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Low with aspirin" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:4]), c("Age \u2265 60", "Australia/New Zealand")))
Sample distribution
pdat <- acs_itt_concurc3_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C3"),
      labels = c("Low", "Intermediate", "Low\nwith aspirin")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
p

Sample distribution
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive-concurrent-stdaspirin.pdf"), p, height = 2.5, width = 6)
Posterior contrast
plot_or_densities(mdrws$compare)

Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-concurrent-stdaspirin-summary-table")
odds_ratio_summary_table_rev(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.32 (0.68, 2.61) 1.40 (0.50) 0.79
Low with aspirin 1.13 (0.60, 2.19) 1.20 (0.41) 0.65
Age ≥ 60 0.51 (0.29, 0.89) 0.53 (0.16) 0.01
Australia/New Zealand 3.15 (0.90, 14.44) 4.28 (3.96) 0.96
Posterior summaries for model parameters.
Code
mfit$summary()
# A tibble: 38 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -3.30e+2 -3.30e+2 3.19    3.10e+0 -3.36e+2 -3.25e+2  1.00    7709.
 2 beta_raw[… -1.05e-1 -1.03e-1 0.247   2.43e-1 -5.18e-1  3.00e-1  1.00   21184.
 3 beta_raw[… -1.65e-1 -1.67e-1 0.235   2.36e-1 -5.53e-1  2.25e-1  1.00   23187.
 4 beta_raw[… -2.70e-1 -2.71e-1 0.115   1.15e-1 -4.57e-1 -8.10e-2  1.00   15091.
 5 beta_raw[…  1.18e+0  1.15e+0 0.710   7.06e-1  8.24e-2  2.42e+0  1.00   19514.
 6 p[1]        3.19e-2  3.14e-2 0.00662 6.52e-3  2.18e-2  4.34e-2  1.00   16106.
 7 p[2]        2.02e-3  1.72e-3 0.00139 1.19e-3  3.96e-4  4.72e-3  1.00   20516.
 8 p[3]        1.08e-3  7.76e-4 0.00103 7.48e-4  7.35e-5  3.15e-3  1.00   17492.
 9 p[4]        1.08e-3  7.71e-4 0.00103 7.55e-4  7.38e-5  3.10e-3  1.00   16262.
10 p[5]        1.08e-3  7.81e-4 0.00101 7.52e-4  7.52e-5  3.10e-3  1.00   15436.
# … with 28 more rows, and 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 0.9718 0.9143 0.9957 0.9689 0.9794 0.9462 1.0284 0.9763
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])

Therapeutic dose

Descriptive summary table
save_tex_table(
  make_summary_table(acs_itt_concurc4_dat, "latex"), 
  file.path("outcomes", "secondary", "7-5-anticoagulation-concurrent-therapeutic-summary"))
make_summary_table(acs_itt_concurc4_dat)
Anticoagulation intervention Patients Known Deaths Any ventilation DAFV, Median (Q1, Q3)
Low-dose 79 75 3 (4%) 6 (8%) 28 (28, 28)
Intermediate-dose 65 63 1 (2%) 3 (5%) 28 (28, 28)
Therapeutic-dose 50 50 6 (12%) 7 (14%) 28 (28, 28)
Overall 194 188 10 (5%) 16 (8%) 28 (28, 28)
Fit model
ctr <- contr.orthonorm
X <- model.matrix( 
  ~ CAssignment + agegte60 + ctry, 
  data = acs_itt_concurc4_nona_dat, 
  contrasts.arg = list(CAssignment = ctr))[, -1]
y_raw <- ordered(acs_itt_concurc4_nona_dat$out_dafv)
y_mod <- as.integer(y_raw)
N <- dim(X)[1]
K <- dim(X)[2]
unique_y <- length(unique(y_mod))
mdat <- list(
  N = N, K = K, J = unique_y, X = X, y = y_mod, y_raw = y_raw,
  p_par = 2 * rep(1 / unique_y, unique_y),
  beta_sd = c(1, 1, 2.5, 1, 1)
)

snk <- capture.output(
  mfit <- ordmod0$sample(
    data = mdat,
    seed = 49135,
    chains = 8,
    parallel_chains = 8,
    iter_warmup = 1000,
    iter_sampling = 2500,
    refresh = 0
))
mdrws <- as_draws_rvars(mfit$draws())
names(mdrws$beta) <- colnames(X)
mdrws$Ccon <- as.vector(ctr(3) %**% mdrws$beta[1:2])
mdrws$Ctrt <- mdrws$Ccon[-1] - mdrws$Ccon[1]
mdrws$compare <- c("Intermediate vs low" = exp(mdrws$Ctrt[1]),
                   "Therapeutic vs low" = exp(mdrws$Ctrt[2]),
                   "Intermediate vs therapeutic" = exp(mdrws$Ctrt[1] - mdrws$Ctrt[2]))
mdrws$OR <- c("Intermediate" = exp(mdrws$Ctrt[1]),
              "Therapeutic" = exp(mdrws$Ctrt[2]),
              setNames(exp(mdrws$beta[3:5]), c("Age \u2265 60", "India", "Australia/New Zealand")))
Sample distribution
pdat <- acs_itt_concurc4_nona_dat %>%
  dplyr::count(
    CAssignment = factor(
      CAssignment, 
      levels = c("C1", "C2", "C4"),
      labels = c("Low", "Intermediate", "Therapeutic")),
    dafv = ordered(out_dafv, levels = 0:28), 
    .drop = F) %>%
  group_by(CAssignment) %>%
  mutate(p = n / sum(n)) %>%
  mutate(cp = cumsum(p)) %>%
  ungroup()
p <- ggplot(pdat, aes(CAssignment, p)) +
  geom_col(aes(fill = fct_rev(dafv))) +
  labs(x = "Anticoagulation intervention", y = "Cumulative proportion") +
  scale_fill_viridis_d("Days alive and\nfree of ventilation", option = "A", direction = -1) +
  guides(fill = guide_legend(reverse = TRUE, ncol = 3)) +
  theme(legend.key.size = unit(0.75, "lines"))
pth <- file.path("outputs", "figures", "outcomes", "secondary")
ggsave(file.path(pth, "outcome-7-5-descriptive-concurrent-therapeutic.pdf"), p, height = 2.5, width = 6)
p

Posterior contrasts
plot_or_densities(mdrws$compare)

Odds ratio summary table
save_tex_table(
  odds_ratio_summary_table_rev(mdrws$OR, "latex"),
  "outcomes/secondary/7-5-primary-model-acs-itt-concurrent-therapeutic-summary-table")
odds_ratio_summary_table_rev(mdrws$OR)
Parameter Median 95% CrI Mean (SD) Pr(OR > 1)
Intermediate 1.46 (0.44, 5.22) 1.81 (1.32) 0.73
Therapeutic 0.57 (0.19, 1.71) 0.67 (0.41) 0.16
Age ≥ 60 0.58 (0.21, 1.61) 0.67 (0.37) 0.15
India 2.63 (0.67, 12.56) 3.64 (3.39) 0.91
Australia/New Zealand 1.13 (0.44, 3.04) 1.28 (0.70) 0.59
Posterior summaries for model parameters.
Code
mfit$summary()
# A tibble: 24 × 10
   variable       mean   median      sd     mad       q5      q95  rhat ess_bulk
   <chr>         <dbl>    <dbl>   <dbl>   <dbl>    <dbl>    <dbl> <dbl>    <dbl>
 1 lp__       -9.86e+1 -9.83e+1 2.47    2.35    -1.03e+2 -95.2     1.00    8006.
 2 beta_raw[… -6.67e-1 -6.60e-1 0.440   0.432   -1.41e+0   0.0417  1.00   20799.
 3 beta_raw[…  6.76e-2  6.60e-2 0.410   0.408   -6.00e-1   0.747   1.00   21069.
 4 beta_raw[… -2.14e-1 -2.15e-1 0.205   0.205   -5.52e-1   0.122   1.00   13273.
 5 beta_raw[…  9.97e-1  9.66e-1 0.751   0.753   -1.96e-1   2.27    1.00   22147.
 6 beta_raw[…  1.26e-1  1.19e-1 0.494   0.489   -6.73e-1   0.946   1.00   16862.
 7 p[1]        4.45e-2  4.13e-2 0.0193  0.0177   1.90e-2   0.0804  1.00   11411.
 8 p[2]        5.75e-3  4.20e-3 0.00543 0.00390  4.74e-4   0.0165  1.00   18867.
 9 p[3]        5.81e-3  4.22e-3 0.00550 0.00397  4.50e-4   0.0165  1.00   15823.
10 p[4]        5.75e-3  4.14e-3 0.00540 0.00390  4.75e-4   0.0164  1.00   15543.
# … with 14 more rows, and 1 more variable: ess_tail <dbl>
Code
mfit$diagnostic_summary()
$num_divergent
[1] 0 0 0 0 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0 0 0 0 0

$ebfmi
[1] 1.0069 1.0749 0.9556 0.9180 0.9832 0.9993 1.0429 0.9973
Code
mcmc_trace(mdrws["beta"])

Code
mcmc_trace(mdrws["alpha"])